Libraries:

# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Functions:

# functions -------------------
source("Functions/veganCovEllipse.R")

Source, tranform, select, and relativize MA & LI veg data (hidden)

What’s different about 5% min plot species for this region?

mali_wider_log <- as.data.frame(lapply(mali_wider, as.logical))

ma.li <- as.data.frame(colSums(mali_wider_log)) %>% 
  rownames_to_column(var = "Species")

# Species from the other 5% calcs
# In NH/ME - no LYQU, CEOR, PAQU, CEAM, ARME, LYBO, RHGL

# For this data set , 5% is ~11 plots:
# Would add: SMGL
# If I stuck to at least 11, would lose: ACRU, COPE, PRSE, PIST, POUN11

rm(mali_wider_log, ma.li)

Source and transform envr & categorical data; merge for predictor df to compare nmds against (hidden)

Run NMDS

# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
                   distance = 'bray',
                   k = 3,
                   trymax = 20,
                   autotransform = FALSE)
## Run 0 stress 0.1395122 
## Run 1 stress 0.1465263 
## Run 2 stress 0.1369332 
## ... New best solution
## ... Procrustes: rmse 0.1183403  max resid 0.3726543 
## Run 3 stress 0.1391444 
## Run 4 stress 0.1458937 
## Run 5 stress 0.1530509 
## Run 6 stress 0.1358958 
## ... New best solution
## ... Procrustes: rmse 0.05522124  max resid 0.1657867 
## Run 7 stress 0.1369326 
## Run 8 stress 0.1358959 
## ... Procrustes: rmse 0.0005646437  max resid 0.001453757 
## ... Similar to previous best
## Run 9 stress 0.1465264 
## Run 10 stress 0.136823 
## Run 11 stress 0.1358959 
## ... Procrustes: rmse 0.0002755069  max resid 0.0008064717 
## ... Similar to previous best
## Run 12 stress 0.1341827 
## ... New best solution
## ... Procrustes: rmse 0.0594534  max resid 0.2139219 
## Run 13 stress 0.1376197 
## Run 14 stress 0.1376205 
## Run 15 stress 0.1391435 
## Run 16 stress 0.1395122 
## Run 17 stress 0.135896 
## Run 18 stress 0.1458936 
## Run 19 stress 0.1369321 
## Run 20 stress 0.1405414 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax
# converges, stress of 0.13 - 0.15

# 2 dimensions
mali_2d <- metaMDS(rel_mali5, 
                   distance = 'bray', 
                   k = 2,
                   try = 40,
                   trymax = 40,
                   autotransform = FALSE)
## Run 0 stress 0.2239971 
## Run 1 stress 0.2187381 
## ... New best solution
## ... Procrustes: rmse 0.1215493  max resid 0.3225023 
## Run 2 stress 0.224739 
## Run 3 stress 0.2047202 
## ... New best solution
## ... Procrustes: rmse 0.1516674  max resid 0.3370896 
## Run 4 stress 0.2066462 
## Run 5 stress 0.2360426 
## Run 6 stress 0.2064071 
## Run 7 stress 0.2216095 
## Run 8 stress 0.2097745 
## Run 9 stress 0.2274517 
## Run 10 stress 0.2191102 
## Run 11 stress 0.2047202 
## ... New best solution
## ... Procrustes: rmse 0.00158856  max resid 0.005415167 
## ... Similar to previous best
## Run 12 stress 0.2150753 
## Run 13 stress 0.2080347 
## Run 14 stress 0.2211264 
## Run 15 stress 0.3659444 
## Run 16 stress 0.2150752 
## Run 17 stress 0.2289216 
## Run 18 stress 0.2272952 
## Run 19 stress 0.2245979 
## Run 20 stress 0.2234159 
## Run 21 stress 0.2132488 
## Run 22 stress 0.224739 
## Run 23 stress 0.2150117 
## Run 24 stress 0.2247403 
## Run 25 stress 0.2232724 
## Run 26 stress 0.2257935 
## Run 27 stress 0.2203333 
## Run 28 stress 0.2097744 
## Run 29 stress 0.2132487 
## Run 30 stress 0.2297719 
## Run 31 stress 0.2185977 
## Run 32 stress 0.2237861 
## Run 33 stress 0.2301115 
## Run 34 stress 0.2261155 
## Run 35 stress 0.2308032 
## Run 36 stress 0.2185277 
## Run 37 stress 0.2150117 
## Run 38 stress 0.2180135 
## Run 39 stress 0.2132487 
## Run 40 stress 0.2146601 
## *** Best solution repeated 1 times
# stress is too high to represent in 2 dimensions

Make correlation data for bi-plots

ma.veg_enr <- envfit(mali_3d, ma.meta.data_veg)
ma.veg_spc.envr <- envfit(mali_3d, rel_mali5)

# site data -------------------
ma.site.scrs <- as.data.frame(scores(mali_3d, display = "sites"))

ma.site.scrs <- cbind(ma.site.scrs, Treat_Type = ma.meta.data_veg$Treat_Type)
ma.site.scrs <- cbind(ma.site.scrs, Region = ma.meta.data_veg$Region)

# species data   -------------------
ma.spp.scrs <- as.data.frame(scores(ma.veg_spc.envr, display = "vectors"))
ma.spp.scrs <- cbind(ma.spp.scrs, Species = rownames(ma.spp.scrs))
ma.spp.scrs <- cbind(ma.spp.scrs, pval = ma.veg_spc.envr$vectors$pvals)

sig.ma_spp.scrs <- subset(ma.spp.scrs, pval<=0.05)

# envr data -------------------
ma.envr.scrs <- as.data.frame(scores(ma.veg_enr, display = "vectors"))
ma.envr.scrs <- cbind(ma.envr.scrs, env.variables = rownames(ma.envr.scrs))
ma.envr.scrs <- cbind(ma.envr.scrs, pval = ma.veg_enr$vectors$pvals)
ma.sig_envr.scrs <- subset(ma.envr.scrs, pval<=0.05)

Need to think about how to represent 3D

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
fort.mali3 <- fortify(mali_3d) %>% 
  filter(score == "sites")



plot_ly(x = fort.mali3$NMDS1, y = fort.mali3$NMDS2, z= fort.mali3$NMDS3, type = "scatter3d", mode = "markers", color = ma.meta.data_veg$Treat_Type, symbols = ma.meta.data_veg$Region)

Run PerMANOVA to look at the impacts of region and treatment on community composition

summary(as.factor(ma.meta.data_veg$Treat_Type))
##  Control   FallRx  Harvest    MowRx SpringRx 
##        5        3        6        3        6
ma.tt <- adonis2(rel_mali5 ~ Treat_Type*Region, data = ma.meta.data_veg)

print(ma.tt)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = rel_mali5 ~ Treat_Type * Region, data = ma.meta.data_veg)
##                   Df SumOfSqs      R2      F Pr(>F)    
## Treat_Type         4   1.7888 0.26616 1.8764  0.002 ** 
## Region             1   0.7921 0.11786 3.3236  0.001 ***
## Treat_Type:Region  2   0.5649 0.08406 1.1852  0.261    
## Residual          15   3.5749 0.53192                  
## Total             22   6.7207 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#each are significant, but the interaction is not